Gene regulation by gonadal hormone receptors underlies brain sex differences

Oestradiol establishes neural sex differences in many vertebrates1–3 and modulates mood, behaviour and energy balance in adulthood4–8. In the canonical pathway, oestradiol exerts its effects through the transcription factor oestrogen receptor-α (ERα)9. Although ERα has been extensively characterized in breast cancer, the neuronal targets of ERα, and their involvement in brain sex differences, remain largely unknown. Here we generate a comprehensive map of genomic ERα-binding sites in a sexually dimorphic neural circuit that mediates social behaviours. We conclude that ERα orchestrates sexual differentiation of the mouse brain through two mechanisms: establishing two male-biased neuron types and activating a sustained male-biased gene expression program. Collectively, our findings reveal that sex differences in gene expression are defined by hormonal activation of neuronal steroid receptors. The molecular targets we identify may underlie the effects of oestradiol on brain development, behaviour and disease.

In mammals, gonadal steroid hormones regulate sex differences in neural activity and behaviour. These hormones establish sex-typical neural circuitry during critical periods of development and activate the display of innate social behaviours in adulthood. Among these hormones, oestradiol is the principal regulator of brain sexual differentiation in mice. In males, the testes briefly activate at birth, generating a sharp rise in testosterone that subsides within hours 10 . Neural aromatase converts circulating testosterone to 17β-oestradiol, which acts through ERα in discrete neuronal populations to specify sex differences in cell number and connectivity 1,3,11 . Despite extensive characterization of the neural circuits controlling sex-typical behaviours 12,13 , the underlying genomic mechanisms by which steroid hormone receptors act in these circuits remain unknown. Recent advancements in low-input and single-cell chromatin profiling methods have provided transformative insights into how transcription factors (TFs) regulate gene expression in small numbers of cells 14 . We set out to use these methods to discover the neuronal genomic targets of ERα and how they coordinate brain sexual differentiation.

Genomic targets of ERα in the brain
To determine the genomic targets of ERα in the brain, we used an established hormone starvation and replacement paradigm that reproducibly elicits sex-typical behaviours 2 and replicates the medium conditions required to detect ERα genomic binding in cell lines 15 . At 4 h after treatment with oestradiol benzoate (E2) or vehicle control, we profiled ERα binding in three interconnected limbic brain regions in which ERα regulates sex-typical behaviours: the posterior bed nucleus of the stria terminalis (BNSTp), medial pre-optic area and posterior medial amygdala 11,12,16 (Fig. 1a). We used the low-input TF profiling method CUT&RUN, which we first validated in MCF-7 breast cancer cells by comparing to a previous dataset for chromatin immunoprecipitation with sequencing (ChIP-seq) of ERα (Extended Data Fig. 1).
We detected 1,930 E2-induced ERα-bound loci in the brain (Fig. 1b, Extended Data Fig. 2 and Supplementary Table 1). The most enriched TF-binding motif in these peaks was the oestrogen response element (ERE), the canonical binding site of oestrogen receptors (Extended Data Fig. 2c, d). Comparison of these ERα-binding sites to those previously detected in peripheral mouse tissues revealed that most are specific to the brain ( Fig. 1c and Extended Data Fig. 2f). Brain-specific ERα binding events were uniquely enriched for synaptic and neurodevelopmental disease Gene Ontology terms, including neurotransmitter receptors, ion channels, neurotrophin receptors and extracellular matrix genes  Table 1). We also found evidence supporting direct crosstalk between oestradiol and neuroprotection, as ERα directly binds loci for the neurotrophin receptors Ntrk2 (also known as Trkb) and Ntrk3 (Extended Data Fig. 2k and Supplementary Table 1). Moreover, ERα targets the genes encoding androgen and progesterone receptors (Ar and Pgr; Supplementary Table 1).
To determine the effects of ERα binding on gene expression and chromatin state, we focused on a single brain region, the BNSTp, given its central role in the regulation of sex-typical behaviours. The BNSTp receives olfactory input through the accessory olfactory bulb and projects to the medial pre-optic area, medial amygdala, hypothalamus and mesolimbic reward pathway 11,17 . We used our oestradiol treatment paradigm and performed translating ribosome affinity purification (TRAP), followed by RNA sequencing (RNA-seq), on the BNSTp from Esr1 Cre/+ ;Rpl22 HA/+ mice, enabling selective capture of ribosome-bound transcripts from Esr1+ cells. We identified 358 genes regulated by oestradiol, including genes known to be induced by E2 in breast cancer, such as Pgr and Nrip1 ( Fig. 1e and Supplementary Table 2). We then validated several of these E2-regulated genes by in situ hybridization (Fig. 1f, Extended Data Fig. 3 and Extended Data Table 1). Genes that contribute to neuron wiring (Brinp2, Unc5b and Enah) and synaptic plasticity (Rcn1 and Irs2) were robustly induced by oestradiol in the Article BNSTp, illustrating how oestradiol signalling may sculpt sexual differentiation of BNSTp circuitry.
To identify oestradiol-responsive chromatin regions, which may involve signalling pathways other than direct ERα binding 18 , we used our oestradiol treatment paradigm and performed assay for transposase-accessible chromatin with sequencing (ATAC-seq) on BNSTp Esr1+ cells collected from Esr1 Cre/+ ;Sun1-GFP lx/+ mice. Across sexes, we detected 7,293 chromatin regions that increase accessibility with treatment (E2-open) as well as 123 regions that decrease accessibility (E2-close; Fig. 1g, Extended Data Fig. 4a-e and Supplementary Table 3). Motif enrichment analysis of these E2-open regions, which occurred primarily at distal enhancer elements (Extended Data Fig. 4c), showed that 89% contain an ERE (Extended Data Fig. 4f), consistent with the observation that nearly all ERα-binding sites overlapped an E2-open region (Extended Data Fig. 4g). These results indicate that direct oestrogen receptor binding, rather than indirect signalling pathways, drives most E2-responsive chromatin regions in the BNSTp 19 .
Although most oestradiol regulation events were shared between sexes in our treatment paradigm, we noted certain sex-dependent effects. Pairwise comparison by sex revealed nearly 300 differential genes between females and males in our TRAP RNA-seq data (Supplementary Table 2). Moreover, we observed 306 genes with a differential response to oestradiol between sexes (Extended Data Fig. 5e, f and Supplementary Table 2). These sex-dependent, E2-responsive genes lacked b 0  Fig. 5g), which may indicate further oestradiol regulation at the translational level 21 . Likewise, across ERα CUT&RUN and ATAC-seq modalities, we observed negligible sex differences and sex-dependent, E2-responsive loci (Extended Data Fig. 5h-j and Supplementary Table 3), demonstrating that females and males mount a similar genomic response to exogenous oestradiol on removal of the hormonal milieu.

Sex differences in gene regulation
Across rodents and humans, the BNSTp of males is approximately 1.5-2 times larger than that of females 22,23 . In mice, this structural dimorphism arises from male-specific neonatal ERα activation, which promotes neuron survival 24,25 . Although BNSTp Esr1+ neurons are known to be GABAergic 16 , the identity of male-biased GABAergic neuron types remains unclear. To characterize these cells, we reanalysed a single-nucleus RNA-seq (snRNA-seq) dataset collected from the BNST of adult, gonadally intact females and males 26 . Seven BNSTp Esr1+ transcriptomic neuron types emerged from this analysis, and two of these marked by Nfix (i1:Nfix) and Esr2 (i3:Esr2) are more abundant in males than in females (Fig. 2a, b and Extended Data Fig. 6a, b). Although a male bias in Esr2/ERβ-labelled cells is known 27 , Nfix expression has not been described previously in the BNSTp. Immunofluorescent staining confirmed that males have twice as many ERα + Nfix + neurons as females ( Fig. 2c and Extended Data Fig. 6c).
To interpret the functional relevance of BNSTp Esr1+ neuron types, we compared their gene expression profiles to the mouse cortical and hippocampal single-cell RNA-seq atlas using MetaNeighbor 28,29 . i1:Nfix neurons uniquely matched the identity of Lamp5+ neurogliaform interneurons 30,31 ( Fig. 2d and Extended Data Fig. 6d, e) and also shared markers (Moxd1 and Cplx3; Extended Data Fig. 6b, f, g) with a male-biased neuron type (i20:Gal/Moxd1) in the sexually dimorphic nucleus of the preoptic area (SDN-POA) that is selectively activated during male-typical mating, inter-male aggression and parenting behaviours 32 . Beyond these two genes, i1:Nfix and i20:Gal/Moxd1 neuron types share a transcriptomic identity, consistent with observed Nfix immunofluorescence in both the BNSTp and SDN-POA ( Fig. 2e and Extended Data Fig. 6h). Together, these results define male-biased neurons in the BNSTp and reveal a common Lamp5+ neurogliaform identity between the BNSTp and SDN-POA.
We next examined sex differences in gene expression and found extensive and robust (false discovery rate < 0.1) sex-biased expression across BNST neuron types (  Table 4). Most sex differences were specific to individual types (for example, Dlg2/PSD-93 and Kctd16 in i1:Nfix neurons), whereas select differences were detected in multiple populations (for example, Tiparp and Socs2; Extended Data Fig. 7b, c). Relative to all other TF genes in the genome, Esr1, along with coexpressed hormone receptors, Ar and progesterone receptor (Pgr), correlated best with sex-biased gene expression ( Fig. 2g and Extended Data Fig. 7e, f), indicating potential regulatory function.
To identify chromatin regions controlling sex differences in BNSTp gene expression, we performed ATAC-seq on BNSTp Esr1+ cells collected from gonadally intact Esr1 Cre/+ ;Sun1-GFP lx/+ mice. Approximately 18,000 regions differed in accessibility between sexes; moreover, these regions localized at sex-biased genes detected in Esr1+ neuron types (Fig. 2h, Extended Data Fig. 7g, h and Supplementary Table 5). By contrast, gonadectomy reduced the number of sex-biased regions to 71 (Fig. 2h and Supplementary Table 5). We compared chromatin accessibility across sexes and gonadal hormone status using k-means clustering and discovered male-specific, but not female-specific, responses to gonadectomy ( Fig. 2i and Extended Data Fig. 7i-k). Notably, chromatin regions that close specifically in males on gonadectomy (cluster 1) primarily contained the androgen response element, whereas regions closing across both sexes (cluster 2) were enriched for the ERE and strongly overlapped E2-open regions (Fig. 2i, j). Thus, in the BNSTp, oestradiol maintains chromatin in an active state across both sexes, whereas testosterone promotes chromatin activation and repression in males. Collectively, these data indicate that gonadal hormone receptors drive adult sex differences in gene expression, largely as a consequence of acute hormonal state.

ERα drives neonatal chromatin state
Sexual dimorphism in BNSTp wiring emerges throughout a 2-week window following birth, well after neural oestradiol has subsided in males. To determine the genomic targets of the neonatal surge, we performed ATAC-seq on BNSTp Esr1+ cells at postnatal day 4 (P4), which corresponds to the onset of male-biased BNSTp cell survival and axonogenesis 33,34 . We detected about 2,000 sex differences in chromatin loci at this time, and nearly all sex differences were dependent on neonatal oestradiol (NE; Fig. 3a, Extended Data Fig. 8a Previous studies have proposed that adult sex differences in behaviour arise from permanent epigenomic modifications induced during the neonatal hormone surge 35 . Our datasets allowed us to examine whether chromatin regions regulated by neonatal hormone maintain sex-biased accessibility into adulthood. Only a small proportion of NE-regulated regions (about 10%) are maintained as sex biased in gonadally intact adults (Extended Data Fig. 9a), implying substantial reprogramming of sex differences as a result of hormonal production during puberty (Fig. 2h). Notably, although most NE-open loci did not maintain male-biased accessibility after puberty, they still localized at adult male-biased genes and clustered around adult male-biased ATAC peaks (Extended Data Fig. 9b-d). These results suggest that certain male-biased genes undergo sequential regulation by ERα and AR in early life and adulthood, respectively.

Sustained sex-biased gene expression
Our identification of approximately 2,000 chromatin regions controlled by the neonatal hormone surge suggests that ERα drives extensive sex differences in the expression of genes that control brain sexual differentiation. To identify these genes, and assess the longevity of their expression, we performed single-nucleus multiome (RNA and ATAC) sequencing on female and male BNST Esr1+ cells collected at P4 and P14, after the closure of the neonatal critical period 3 ( Fig. 3b and Extended Data Fig. 10a, b). We profiled 14,836 cells and found that Esr1+ neuron identity is largely the same across P4, P14 and adulthood 36 ( Fig. 3b and Extended Data Fig. 10c-f).
To identify TFs regulating Esr1+ neuron identity, we ranked TFs on their potential to control chromatin accessibility and their expression specificity across neuron types 37 (Extended Data Fig. 10g). This approach uncovered canonical GABAergic identity TF genes a priori, including Lhx6, Prox1 and Nkx2-1, as well as regulators Zfhx3 and Nr4a2 (Extended Data Fig. 10g). In addition, Nfix was predicted to regulate the identity of the male-biased i1:Nfix neuron type ( Fig. 3c and Extended Data Fig. 10f, g). Profiling Nfix binding in the adult BNSTp confirmed that the binding sites of this factor, including at the Nfix locus itself, are maintained in an active state primarily in i1:Nfix neurons (Fig. 3c (Fig. 3d). These data suggest that, in addition to specifying the chromatin landscape, identity TFs may dictate the cellular response to neonatal oestradiol by influencing ERα binding.
Differential expression analysis across Esr1+ neuron types on P4 identified >400 sex-biased genes (Fig. 3e, Extended Data Fig. 11a and Supplementary Table 9). Performing RNA-seq on BNSTp Esr1+ cells collected from females treated at birth with vehicle or oestradiol showed that these sex differences largely arise as a consequence of the neonatal surge (Extended Data Fig. 11b and chromatin state occurred in neurons lacking Cyp19a1/aromatase expression ( Fig. 3e-g), indicative of non-cell-autonomous oestradiol signalling.
To link our chromatin and gene expression data, we constructed a gene regulatory map across Esr1+ neuron types consisting of sex-biased genes and NE-regulated enhancers with correlated accessibility (Fig. 3e and Extended Data Fig. 11f, g). This map demonstrates both divergent responses across neuron types as well as neuron-type-specific enhancers for common sex-biased targets. Notably, we identified Arid1b, an autism spectrum disorder candidate gene, among genes regulated by distinct enhancers across neuron types (Extended Data Fig. 11g).
Further examination showed that about 40% of high-confidence (family-wise error rate ≤ 0.05) autism spectrum disorder candidate genes 38 , including Grin2b, Scn2a1 (also known as Scn2a) and Slc6a1, contained NE-open chromatin regions and ERα occupancy (Extended Data Fig. 8j and Supplementary Table 6).
We also examined whether sex-biased genes, and their corresponding enhancers, are sustained across the neonatal critical period by comparing Esr1+ neurons between P4 and P14. Although the total number of sex-biased genes declined between P4 and P14, a subset persisted as sex biased throughout the neonatal critical window ( Fig. 3h Lrp1b 660 kb  Article expressed genes on P4 persisted as sex biased on P14. These genes regulate distinct components of neural circuit development, including neurite extension (Klhl1 and Pak7 (also known as Pak5)), axon pathfinding (Epha3 and Nell2), neurotransmission (Kcnab1 and Scg2) and synapse formation (Il1rap and Tenm2; Fig. 3h and Extended Data Fig. 11h). Together, these results show that neonatal ERα activation drives the epigenetic maintenance of a gene expression program that facilitates sexual differentiation of neuronal circuitry.

Sustained sex differences require ERα
The adult display of male mating and territoriality behaviours requires ERα expression in GABAergic neurons 16 . To determine whether ERα is also required for sustained sex differences in gene expression, we performed snRNA-seq on 38,962 BNST GABAergic neurons isolated from P14 conditional mutant males lacking ERα (Vgat Cre ;Esr1 lx/lx ;Sun1-GFP lx ), and littermate control females and males (Vgat Cre ;Esr1 +/+ ;Sun1-GFP lx ; Fig. 4a and Extended Data Fig. 12a-d). GABAergic neurons in ERα-mutant males did not deviate from P14 control or adult BNST neuron types ( Fig. 4a and Extended Data Fig. 12b), indicating that ERα is dispensable for neuron identity. However, the abundance of male-biased i1:Nfix and i3:Esr2 neurons dropped to female levels in Vgat Cre ;Esr1 lx/lx males ( Fig. 4b and Extended Data Fig. 12d), suggesting that neonatal ERα activation is essential for their male-typical abundance. Differential expression analysis between control females and control or conditional ERα-knockout (KO) males in each neuron type established that ERα is required for nearly all sexually dimorphic gene expression, with the exception of those located on the Y chromosome or escaping X inactivation (Fig. 4b, Extended Data Fig. 12e and Supplementary Table 11). Notably, ERα-KO males exhibited feminized expression of sex-biased genes ( Fig. 4c and Extended Data Fig. 12f). Together, these findings demonstrate that the neonatal hormone surge drives a sustained male-typical gene expression program through activation of a master regulator TF, ERα (Fig. 4c).

Discussion
Here we identify the genomic targets of ERα in the brain and demonstrate that BNSTp sexual differentiation is defined by both male-biased cell number and gene expression. We find that sexual dimorphism in the BNSTp equates to increased numbers of i1:Nfix and i3:Esr2 neurons in males. The transcriptomic identity of i1:Nfix neurons resembles that of cortical Lamp5+ neurogliaform interneurons, which provide regional inhibition through synaptic and ambient release of GABA 39 . As all BNSTp neurons and much of the POA are GABAergic, we predict that higher numbers of i1:Nfix inhibitory neurons enables stronger disinhibition of downstream projection sites. The net result is a gain of responses to social information, leading to male-typical levels of mounting or attacking 40 . Male-biased populations of inhibitory neurons also modulate sex-typical behaviours in Drosophila, but they do not rely on gonadal hormones to specify sex-biased enhancers 41-45 . In vertebrates, hormone receptor signalling may have evolved to coordinate gene regulation throughout a neural circuit as a strategy for controlling context-dependent behavioural states. Moreover, the association between hormone receptor target genes identified here and human neurological and neurodevelopmental conditions may explain the notable sex biases of these diseases.
Our data show that the neonatal hormone surge activates ERα to drive a sustained male-biased gene expression program in the developing brain. We speculate that this program establishes male-typical neuronal connectivity across the neonatal critical period and potentially primes the response to hormone receptor activation at puberty. In the adult brain, gonadectomy ablated sex differences in chromatin accessibility, and under these conditions, Esr1+ neurons of both sexes exhibited a similar genomic response to exogenous oestradiol. Together, these findings suggest that although sex differences in developmental gonadal hormone signalling establish dimorphisms in BNSTp circuitry, the genome remains responsive to later alterations in the hormonal milieu. Likewise, manipulating hormonal status, circuit function or individual genes consistently demonstrates that both sexes retain the potential to engage in behaviours typical of the opposite sex 46-49 . This study implicates puberty as a further critical period for sexual differentiation of gene regulation and provides an archetype for studying hormone receptor action across life stages, brain regions and species.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-04686-1.

Animals
All animals were maintained on a 12-h light/12-h dark cycle and provided food and water ad libitum. All mouse experiments were performed under strict guidelines set forth by the CSHL Institutional Animal Care and Use Committee. All animals were randomly assigned to experimental groups. Esr1 Cre (ref. 50 ), Rpl22 HA (ref. 51 ), ROSA26 CAG-Sun1-sfGFP-Myc (ref. 52 ); abbreviated as Sun1-GFP), Vgat Cre (ref. 53 ) and C57Bl6/J wild-type mice were obtained from Jackson Labs. Esr1 lx mice were received from S. A. Khan 54 . Adult male and female mice were used between 8 and 12 weeks of age. For adult hormone treatment experiments, animals were euthanized for tissue collection 4 h after subcutaneous administration of 5 µg E2 (Sigma E8515) suspended in corn oil (Sigma C8267) or vehicle 3 weeks post-gonadectomy. For neonatal CUT&RUN, ATAC-seq and RNA-seq experiments, animals were treated with 5 µg E2 or vehicle on P0 and collected 4 h later (ERα CUT&RUN) or 4 days later (ATAC-seq and nuclear RNA-seq). For neonatal multiome, snRNA-seq and IF quantification, animals were collected on P4 (multiome) or P14 (multiome, snRNA-seq and IF staining).

Cell lines
Cell lines include mHypoA clu-175 clone (Cedarlane Labs) and MCF-7 (ATCC). Cell lines were not tested for mycoplasma contamination. Cells were maintained in standard DMEM supplemented with 10% FBS and penicillin/streptomycin. Before CUT&RUN, MCF7 cells were grown in phenol-red-free DMEM medium containing 10% charcoal-stripped FBS and penicillin/streptomycin for 48 h and then treated with 20 nM 17-β-oestradiol or vehicle (0.002% ethanol) for 45 min.

Adult RNA-seq and in situ hybridization
Experiments were performed as previously described 55 . Briefly, the BNSTp was microdissected following rapid decapitation of deeply anaesthetized adult Esr1 Cre/+ ;Rpl22 HA/+ mice. Tissue homogenization, immunoprecipitation and RNA extraction were performed, and libraries were prepared from four biological replicate samples (each consisting of 8-9 pooled animals) using NuGEN Ovation RNA-Seq kits (7102 and 0344). Multiplexed libraries were sequenced with 76-bp single-end reads on the Illumina NextSeq. Validation by in situ hybridization staining and quantification was performed by an investigator blinded to experimental condition, as previously described 16,55 . Riboprobe sequences are listed in Extended Data Table 1.

Isolation of nuclei from adult mice for ATAC-seq
Adult Esr1 Cre/+ ;Sun1-GFP lx/+ mice (four pooled per condition) were deeply anaesthetized with ketamine/dexmedetomidine. Sections of 500 µm spanning the BNSTp were collected in an adult mouse brain matrix (Kent Scientific) on ice. The BNSTp was microdissected and collected in 1 ml of cold supplemented homogenization buffer (250 mM sucrose, 25 mM KCl, 5 mM MgCl 2 , 120 mM tricine-KOH, pH 7.8), containing 1 mM dithiothreitol, 0.15 mM spermine, 0.5 mM spermidine and 1× EDTA-free PIC (Sigma Aldrich 11873580001). The tissue was dounce homogenized 15 times in a 1-ml glass tissue grinder (Wheaton) with a loose pestle. Next, 0.3% IGEPAL CA-630 was added, and the suspension was homogenized five times with a tight pestle. The homogenate was filtered through a 40-µm strainer and then centrifuged at 500g for 15 min at 4 °C. The pellet was resuspended in 0.5 ml homogenization buffer containing 1 mM dithiothreitol, 0.15 mM spermine, 0.5 mM spermidine and 1× EDTA-free PIC. A total of 30,000 GFP + nuclei were collected into cold ATAC-RSB (10 mM Tris-HCl pH 7.5, 10 mM NaCl, 3 mM MgCl 2 ) using the Sony SH800S Cell Sorter (purity mode) with a 100-µm sorting chip. After sorting, 0.1% Tween-20 was added, and the nuclei were centrifuged at 500g for 5 min at 4 °C. The pellet of nuclei was directly resuspended in transposition reaction mix.

ATAC-seq library preparation
Tn5 transposition was performed using the OMNI-ATAC protocol 56 . A 2.5 µl volume of Tn5 enzyme (Illumina 20034197) was used in the transposition reaction. Libraries were prepared with NEBNext High-Fidelity 2× PCR Master Mix (NEB M0541L), following the standard protocol. After the initial five cycles of amplification, another four cycles were added, on the basis of qPCR optimization. Following amplification, libraries were size selected (0.5×-1.8×) twice with AMPure XP beads (Beckman Coulter A63880) to remove residual primers and large genomic DNA. Individually barcoded libraries were multiplexed and sequenced with paired-end 76-bp reads on an Illumina NextSeq, using either the Mid or High Output Kit.

Cell line CUT&RUN
To collect cells for CUT&RUN, cells were washed twice with Hank's buffered salt solution (HBSS) and incubated for 5 min with pre-warmed 0.5% trypsin-EDTA (10×) at 37 °C/5% CO 2 . Trypsin was inactivated with DMEM supplemented with 10% FBS and penicillin/streptomycin (mHypoA cells) or phenol-red-free DMEM supplemented with 10% charcoal-stripped FBS and penicillin/streptomycin (MCF-7 cells). After trypsinizing, cells were centrifuged at 500g in a 15-ml conical tube and resuspended in fresh medium. CUT&RUN was performed as previously described 14 , with minor modifications. Cells were washed twice in wash buffer (20 mM HEPES, pH 7.5, 150 mM NaCl, 0.5 mM spermidine, 1× PIC, 0.02% digitonin). Cell concentration was measured on a Countess II FL Automated Cell Counter (Thermo Fisher). A total of 25,000 cells were used per sample. Cells were bound to 20 µl concanavalin A beads (Bangs Laboratories, BP531), washed twice in wash buffer, and incubated overnight with primary antibody (ERα: Santa Cruz sc-8002 or EMD Millipore Sigma 06-935, Nfix: Abcam ab101341) diluted 1:100 in antibody buffer (wash buffer containing 2 mM EDTA). The following day, cells were washed twice in wash buffer, and 700 ng ml −1 protein A-MNase (pA-MNase, prepared in-house) was added. After 1 h incubation at 4 °C, cells were washed twice in wash buffer and placed in a metal heat block on ice. pA-MNase digestion was initiated with 2 mM CaCl 2 . After 90 min, digestion was stopped by mixing 1:1 with 2× stop buffer (340 mM NaCl, 20 mM EDTA, 4 mM EGTA, 50 µg ml −1 RNase A, 50 µg ml −1 glycogen, 0.02% digitonin). Digested fragments were released by incubating at 37 °C for 10 min, followed by centrifuging at 16,000g for 5 min at 4 °C. DNA was purified from the supernatant by phenol-chloroform extraction, as previously described 14 .

CUT&RUN library preparation
Cell line CUT&RUN libraries were prepared using the SMARTer Thru-PLEX DNA-seq Kit (Takara Bio R400676), with the following PCR conditions: 72 °C for 3 min, 85 °C for 2 min, 98 °C for 2 min, (98 °C for 20 s, 67 °C for 20 s, 72 °C for 30 s) × 4 cycles, (98 °C for 20 s, 72 °C for 15 s) × 14 cycles (MCF7) or 10 cycles (mHypoA). Brain CUT&RUN libraries were prepared using the same kit with 10 PCR cycles. All samples were size selected with AMPure XP beads (0.5×-1.7×) to remove residual adapters and large genomic DNA. Individually barcoded libraries were multiplexed and sequenced with paired-end 76-bp reads on an Illumina NextSeq, using either the Mid or High Output Kit. For the mHypoA experiment, samples were sequenced with paired-end 25-bp reads on an Illumina MiSeq.

Neonatal bulk ATAC-seq
Female and male Esr1 Cre/+ ;Sun1-GFP lx/+ mice were injected subcutaneously with 5 µg E2 or vehicle on P0 and collected on P4 (4-5 animals pooled per condition and per replicate). The BNSTp was microdissected, as described above, and collected in 300 µl of cold, supplemented homogenization buffer. Nuclei were extracted as described for the adult brain. After filtering through a 40-µm strainer, the nuclei were diluted 3:1 with 600 µl of cold, supplemented homogenization buffer and immediately used for sorting. A total of 30,000 GFP + nuclei were collected into cold ATAC-RSB buffer using the Sony SH800S Cell Sorter (purity mode) with a 100-µm sorting chip. After sorting, nuclei transposition and library preparation were performed, as described above.

P0 ERα CUT&RUN
Female Esr1 Cre/+ ;Sun1-GFP lx/+ mice were injected subcutaneously with 5 µg E2 or vehicle on P0 and collected 4 h later (5 animals pooled per condition and per replicate). The BNSTp, MPOA and MeA were microdissected, and nuclei were extracted, as described for the neonatal bulk ATAC-seq experiment. After filtering through a 40-µm strainer, the nuclei were diluted 3:1 with 600 µl of cold, supplemented homogenization buffer. A 2 mM concentration of EDTA was added, and the sample was immediately used for sorting. A total of 150,000 GFP + nuclei were collected into cold CUT&RUN wash buffer using the Sony SH800S Cell Sorter (purity mode) with a 100-µm sorting chip. GFP − events were collected into cold CUT&RUN wash buffer, and 150,000 nuclei were subsequently counted on the Countess II FL Automated Cell Counter for ERα− and IgG negative-control CUT&RUN. All subsequent steps were performed as described for the adult brain CUT&RUN experiments. P0 CUT&RUN libraries were prepared with 10 PCR cycles.

Neonatal single-nucleus multiome sequencing
The BNST was microdissected fresh from P4 and P14 female and male Esr1 Cre/+ ;Sun1-GFP lx/+ mice, as described above (4-5 animals pooled per condition). Nuclei were extracted and prepared for sorting, as performed for the neonatal bulk ATAC-seq experiment, with the inclusion of 1 U µl −1 Protector RNase inhibitor (Sigma) in the homogenization buffer. A total of 40,000-50,000 GFP + nuclei were collected into 1 ml of cold ATAC-RSB buffer, supplemented with 0.1% Tween-20, 0.01% digitonin, 2% sterile-filtered BSA (Sigma A9576) and 1 U µl −1 Protector RNase inhibitor. The nuclei were centrifuged in a swinging-bucket rotor at 500g for 10 min at 4 °C. About 950 µl of supernatant was carefully removed, and 200 µl 10x Genomics dilute nuclei buffer was added to the side of the tube without disturbing the pellet. The nuclei were centrifuged again at 500g for 10 min at 4 °C. About 240 µl of supernatant was carefully removed, and the nuclei were resuspended in the remaining volume (about 7 µl). Samples were immediately used for the 10x Genomics Single Cell Multiome ATAC + Gene Expression kit (1000285), following the manufacturer's instructions. snRNA-seq and snATAC-seq libraries were sequenced on an Illumina NextSeq, using the High Output kit. Each sample was sequenced to a depth of about 40,000-80,000 mean reads per cell for the snATAC library and about 40,000-50,000 mean reads per cell for the snRNA library.

P14 snRNA-seq
The BNSTp was microdissected from P14 female and male Vgat Cre ;Esr1 +/+ ;Sun1-GFP lx and male Vgat Cre ;Esr1 lx/lx ;Sun1-GFP lx mice. Tissue samples from individual animals were immediately flash frozen in an ethanol dry-ice bath and stored at −80 °C until n = 3 animals were collected per group. On the day of the experiment, tissue samples were removed from −80 °C and maintained on dry ice. With the tissue still frozen, cold, supplemented homogenization buffer was added to the tube, and the tissue was immediately transferred to a glass homogenizer and mechanically dounced and filtered, as described for our other neonatal experiments. A total of 80,000-90,000 GFP + nuclei were collected into 100 µl of cold ATAC-RSB buffer, supplemented with 1% sterile-filtered BSA (Sigma A9576), and 1 U µl −1 Protector RNase inhibitor, in a 0.5-ml DNA lo-bind tube (Eppendorf) pre-coated with 30% BSA. After collection, nuclei were pelleted with two rounds of gentle centrifugation (200g for 1 min) in a swinging-bucket centrifuge at 4 °C. After the second round, the supernatant was carefully removed, leaving about 40 µl in the tube. The nuclei were gently resuspended in this remaining volume and immediately used for the 10x Genomics Single Cell 3′ Gene Expression kit v3 (1000424), following the manufacturer's instructions. Each biological sample was split into two 10× lanes, producing 6 libraries that were pooled and sequenced on an Illumina NextSeq 2000 to a depth of about 45,000-60,000 mean reads per cell.

Neonatal nuclear RNA-seq
Female Esr1 Cre/+ ;Sun1-GFP lx/+ mice were injected subcutaneously with 5 µg E2 or vehicle on P0. Four days later, animals were rapidly decapitated, and 400-µm sections were collected in cold homogenization buffer using a microtome (Thermo Scientific Microm HM 650V). The BNST was microdissected (4 animals pooled per condition) and collected in 1 ml of cold, supplemented homogenization buffer containing 0.4 U ml −1 RNAseOUT (Thermo Fisher, 10777019). Nuclei were isolated as described for neonatal bulk ATAC-seq. A total of 12,000 GFP + nuclei were collected into cold Buffer RLT Plus supplemented 1:100 with β-mercaptoethanol (Qiagen, 74034) using the Sony SH800S Cell Sorter (purity mode) with a 100-µm sorting chip. Nuclei lysates were stored at −80 °C until all replicates were collected. Nuclei samples for all replicates were thawed on ice, and RNA was isolated using the Qiagen RNeasy Plus Micro Kit (74034). Strand-specific RNA-seq libraries were prepared using the Ovation SoLo RNA-seq system (Tecan Genomics, 0501-32), following the manufacturer's guidelines. Individually barcoded libraries were multiplexed and sequenced with single-end 76-bp reads on an Illumina NextSeq, using the Mid Output Kit.
CUT&RUN data analysis. CUT&RUN differential peak calling was performed with DiffBind v2.10.0(ref. 64 ). A count matrix was created from individual replicate BAM and MACS2 narrowpeak files (n = 2 per condition). Consensus peaks were recentred to ±100 bp around the point of highest read density (summits = 100). Contrasts between sex and treatment were established (categories = c(DBA_TREATMENT, DBA_CONDITION)), and edgeR 65 was used for differential peak calling. Differential ERα peaks with P adj < 0.1 were used for downstream analysis. For Nfix, differential peaks with a P adj < 0.1 and abs(log 2 [FC]) > 1 were used for downstream analysis. Differential peak calling for the MCF-7 CUT&RUN experiment was performed with DESeq2 (P adj < 0.1) in DiffBind. Differential peak calling for the P0 ERα CUT&RUN experiment was performed with DESeq2 (P adj < 0.01) in DiffBind. To identify sex-dependent, oestradiol-responsive peaks for adult brain ERα CUT&RUN, the DiffBind consensus peakset count matrix was used as input to edgeR, and an interaction between sex and treatment was tested with glmQLFTest.
MCF7 ERα CUT&RUN data were compared to MCF7 ERα ChIP-seq data from ref. 79 (GEO: GSE59530). Single-end ChIP-seq fastq files for two vehicle-treated and two 17β-oestradiol (E2)-treated IP and input samples were accessed from the Sequence Read Archive and processed identically to ERα CUT&RUN data, with the exception of fragment size filtering. Differential ERα ChIP-seq peak calling was performed using DiffBind DESeq2 (P adj < 0.01). DeepTools was used to plot CPM-normalized ERα CUT&RUN signal at E2-induced ERα ChIP-seq binding sites. DREME 80 and AME were used to compare de novo and enriched motifs between E2-induced MCF7 ERα CUT&RUN and ChIP-seq peaks.
Adult RNA-seq data processing and analysis. Reads were adapter trimmed and quality filtered (q > 30) (http://hannonlab.cshl.edu/ fastx_toolkit/), and then mapped to the mm10 reference genome using STAR 81 . The number of reads mapping to the exons of each gene was counted with featureCounts 82 , using the NCBI RefSeq mm10 gene annotation. Differential gene expression analysis was performed using DESeq2 (ref. 83 ) with the following designs: effect of treatment (design = ~ batch + hormone), effect of sex (design = ~ batch + sex), two-way comparison of treatment and sex (design = ~ batch + hormone_sex), four-way comparison (design = ~ 0 + hormone_sex) and sex-treatment interaction (design = ~ batch + sex + hormone + sex:hormone).
Adult GDX treatment ATAC-seq data analysis. ATAC-seq differential peak calling was performed with DiffBind v2.10.0. A DiffBind dba object was created from individual replicate BAM and MACS2 narrowPeak files (n = 3 per condition). A count matrix was created with dba.count, and consensus peaks were recentred to ±250 bp around the point of highest read density (summits = 250). Contrasts between sex and treatment were established (categories = c(DBA_TREATMENT, DBA_CONDITION)), and edgeR was used for differential peak calling. Differential peaks with an FDR < 0.05 and abs(log 2 [FC]) > 1 or abs(log 2 [FC] )> 0 were used for downstream analysis. DeepTools computeMatrix and plotHeatmap were used to plot mean ATAC CPM at E2-open ATAC peaks. To identify sex-dependent, oestradiol-responsive peaks, the DiffBind consensus peakset count matrix was used as input to edgeR, and an interaction between sex and treatment was tested with glmQLFTest. E2-open ATAC peaks and total vehicle or E2 ATAC peaks (intersected across replicate and sex for each treatment condition) were annotated to NCBI RefSeq mm10 genes using ChIPseeker. ClusterProfiler was used to calculate the enrichment of GO biological process terms. DO and HGNC gene family enrichment was performed on E2-open ATAC peak-associated genes, as described above for ERα CUT&RUN analysis. BETA (basic mode, -d 500000) 68 was used to determine whether E2-open ATAC peaks were significantly overrepresented at E2-regulated RNA-seq genes (P < 0.01), as well as sex-dependent E2-regulated genes (P < 0.01), compared to non-differential, expressed genes. Motif enrichment analysis of E2-open ATAC peaks was performed with AME, using the 2020 JASPAR core non-redundant vertebrate database. FIMO was used to determine the percentage of E2-open ATAC peaks containing the enriched motifs shown in Extended Data Fig. 4h, i.
Adult gonadally intact ATAC-seq analysis. ATAC-seq differential peak calling and comparison between gonadally intact (abbreviated as intact) and GDX ATAC samples were performed with DiffBind v2.10.0 and edgeR. A DiffBind dba object was created from individual replicate BAM and MACS2 narrowPeak files for the four groups: female intact (n = 2), male intact (n = 2), female GDX vehicle treated (n = 3), male GDX vehicle treated (n = 3). A count matrix was created with dba.count, and consensus peaks were recentred to ±250 bp around the point of highest read density (summits = 250). The consensus peakset count matrix was subsequently used as input to edgeR. Differential peaks (abs(log 2 [FC]) > 1, P adj < 0.05) were calculated between female intact and male intact and between female GDX vehicle treated and male GDX vehicle -treated groups using glmQLFTest. BETA was used to assess statistical association between gonadally intact, sex-biased ATAC peaks and sex DEGs called in BNSTp Esr1+ snRNA-seq clusters (top 500 genes per cluster, ranked by P adj ). Sex DEGs ranked by ATAC regulatory potential score 68 , a metric that reflects the number of sex-biased peaks and distance of sex-biased peaks to the TSS, are shown in Extended Data Fig. 7g. HGNC gene family enrichment was performed on sex DEGs, using a background of expressed genes in any of the seven BNSTp Esr1+ clusters.
To identify differential peaks across the four conditions, an ANOVA-like design was created in edgeR by specifying multiple coefficients in glmQLFTest (coefficient = 2:4). A matrix of normalized counts in these differential peaks (P adj < 0.01) was clustered using k-means clustering (kmeans function in R), with k = 4 and iter. max = 50. For each k-means cluster, the cluster centroid was computed, and outlier peaks in each cluster were excluded on the basis of having low Pearson's correlation with the cluster centroid (R < 0.8). Depth-normalized ATAC CPM values in these peak clusters are shown in Fig. 2i (mean across biological replicates per group) and Extended Data Fig. 7 (individual biological replicates). Peak cluster overlap with E2-open ATAC loci (abs(log 2 [FC]) > 0, P adj < 0.05) was computed with bedtools intersect (-wa). For each peak cluster, motif enrichment analysis was performed by first generating a background peak list (matching in GC content and accessibility) from the consensus ATAC peak matrix using chromVAR (addGCBias, getBackgroundPeaks) 84 , and then calculating enrichment with AME using the background peak list as the control (--control background peaks). In Fig. 2i, the JASPAR 2020 AR motif (MA0007.3) is labelled as ARE, and the ESR2 motif (MA0258.2) is labelled as ERE.
Adult snRNA-seq and single-cell RNA-seq analysis. Mouse BNST snRNA-seq data containing 76,693 neurons across 7 adult female and 8 adult male biological replicates 26 were accessed from GEO: GSE126836 and loaded into a Seurat object 85 . Mouse MPOA single-cell RNA-seq data containing 31,299 cells across 3 adult female and 3 adult male biological replicates 32 were accessed from GEO: GSE113576 and loaded into a Seurat object. Cluster identity, replicate and sex were added as metadata features to each Seurat object. Pseudo-bulk RNA-seq analysis was performed to identify sex differences in gene expression in the BNST snRNA-seq dataset. Briefly, the Seurat object was converted to a SingleCellExperiment object (as.SingleCellExperiment). Genes were filtered by expression (genes with >1 count in ≥5 nuclei). NCBI-predicted genes were removed. For each cluster, nuclei annotated to the cluster were subsetted from the main Seurat object. Biological replicates containing ≤20 nuclei in the subsetted cluster were excluded. Gene counts were summed for each biological replicate in each cluster. Differential gene expression analysis across sex in each cluster was performed on the filtered, aggregated count matrix using DESeq2 (design = ~ sex) with alpha = 0.1. The BNSTp_Cplx3 cluster was excluded, as none of the replicates in this cluster contained more than 20 nuclei. Clusters containing ≥25% nuclei with ≥1 Esr1 counts in the main Seurat object were classified as Esr1+ (i1:Nfix, i2:Tac2, i3:Esr2, i4:Bnc2, i5:Haus4, i6:Epsti1, i7:Nxph2, i8:Zeb2, i9:Th, i10:Synpo2, i11:C1ql3, i12:Esr1, i13:Avp, i14:Gli3). To identify TFs that correlate with sex DEG number per cluster (Fig. 2g), a linear regression model with percentage of TF expression as the predictor variable and sex DEG number per cluster as the response variable was generated using the lm function in R stats (formula = percentage of TF expression ~ DEG number). This model was tested for all TFs in the SCENIC 86 mm10 database. All TFs were then ranked by R 2 to identify those most predictive of sex DEG number, and the ranked R 2 values are shown in Fig. 2g.
To visualize BNSTp Esr1+ snRNA-seq data (Fig. 2a), BNSTp Esr1+ clusters were subsetted from the main Seurat object. Gene counts were normalized and log transformed (LogNormalize), and the top 2,000 variable features were identified using FindVariableFeatures (selection.method = vst). Gene counts were scaled, and linear dimensionality reduction was performed by principal component analysis (runPCA, npcs = 10). BNSTp Esr1+ clusters were visualized with UMAP (runUMAP, dims = 10). To generate the heatmaps in Extended Data Fig. 7a, pseudo-bulk counts for each biological replicate included in the analysis were normalized and transformed with variance-stabilizing transformation (DESeq2 vst), subsetted for sex-biased genes in each cluster, and z-scaled across pseudo-bulk replicates.
To examine differential abundance of BNSTp Esr1+ clusters between sexes (Fig. 2b), the proportion of total nuclei in each BNSTp Esr1+ cluster was calculated for each biological replicate. After calculating the proportions of nuclei, sample MALE6 was excluded as an outlier for having no detection (0 nuclei) of i1:Nfix and i2:Tac2 clusters and overrepresentation of the i5:Haus4 cluster. The one-sided Wilcoxon rank-sum test (wilcox.test in R stats) was used to test for male-biased abundance of nuclei across biological replicates in each cluster. P values were adjusted for multiple hypothesis testing using the Benjamini-Hochberg procedure (method = fdr).
To identify the enrichment of Lamp5+ subclass markers in BNSTp and MPOA Esr1+ clusters (Extended Data Fig. 6e), a Seurat object was created from the Allen Brain Atlas Cell Types dataset. Gene counts per cell were normalized and log transformed (LogNormalize), and subclass-level marker genes were calculated with the Wilcoxon rank-sum test (Find-AllMarkers, test.use = wilcox, min.diff.pct = 0.2). The mean expression of Lamp5+ subclass markers (avg_log[FC] > 0.75, P adj < 0.05, <40% in non-Lamp5+ subclasses) was calculated in BNSTp and MPOA Esr1+ clusters and visualized using pheatmap.
MetaNeighbor analysis. MetaNeighbor 28 was used to quantify the degree of similarity between BNSTp Esr1+ clusters and MPOA Esr1+ clusters and between BNSTp Esr1+ clusters and cortical/hippocampal GABAergic neuron subclasses from the Allen Brain Atlas Cell Types database 29 . Briefly, the BNST and MPOA Seurat objects were subsetted for Esr1+ clusters, and then transformed and merged into one Single-CellExperiment object. For the BNSTp and cortex comparison, BNSTp Esr1+ clusters were merged into a SingleCellExperiment with cortical/ hippocampal GABAergic cortical clusters. Unsupervised MetaNeighbor analysis was performed between BNST and MPOA clusters, and between BNST and cortical/hippocampal clusters, using highly variable genes identified across datasets (called with the variableGenes function). The median AUROC value per cortical/hippocampal GABAergic subclass across Allen Brain Atlas datasets for each BNSTp Esr1+ cluster is shown in Fig. 2d.
Neonatal bulk ATAC-seq analysis. Differential peak calling on the neonatal bulk ATAC-seq experiment was performed with DiffBind v2.10.0 and edgeR. A count matrix was created from individual replicate BAM and MACS2 narrowpeak files (n = 3 per condition). Consensus peaks were recentred to ±250 bp around the point of highest read density (summits = 250), and the consensus peakset count matrix was subsequently used as input to edgeR. Differential peaks across the three treatment groups (NV female, NV male, NE female) were calculated by specifying multiple coefficients in glmQLFTest (coefficient = 4:5). To identify accessibility patterns across differential peaks (P adj < 0.05), a matrix of normalized counts in differential peaks was hierarchically clustered using pheatmap, and the resulting dendrogram tree was cut with k = 6 to achieve 6 peak clusters (Extended Data Fig. 8a). The two largest clusters were identified as having higher accessibility in NV males and NE females compared to NV females (cluster 3, labelled as NE open), or lower accessibility in NV male and NE female compared to NV females (cluster 5, labelled as NE close). Motif enrichment analysis of NE-open peaks was performed with AME using the 2020 JASPAR core non-redundant vertebrate database. GO biological process, DO and HGNC gene family enrichment analyses were performed, as described above for adult GDX treatment ATAC-seq data analysis.
Neonatal single-nucleus multiome data processing and analysis. Raw sequencing data were processed using the Cell Ranger ARC pipeline (v2.0.0) with the cellranger-arc mm10 reference. Default parameters were used to align reads, count unique fragments or transcripts, and filter high-quality nuclei. Individual HDF5 files for each sample containing RNA counts and ATAC fragments per cell barcode were loaded into Seurat (Read10X_h5). Nuclei with lower-end ATAC and RNA QC metrics (<1,000 ATAC fragments, <500 counts, nucleosomal signal > 3, TSS enrichment < 2) were removed. DoubletFinder 87 was then used to remove predicted doublets from each sample (nExp = 9% of nuclei per sample). Following doublet removal, nuclei surpassing upper-end ATAC and RNA QC metrics (>60,000 ATAC fragments, >20,000 RNA counts, >6,000 genes detected) were removed. After filtering, Seurat objects for each sample were subsetted for the RNA assay and merged. Gene counts were normalized and log transformed (LogNormalize), and the top 2,000 variable features were identified using FindVariableFeatures (selection.method = 'vst'). Gene counts were scaled, regressing out the following variables: number of RNA counts, number of RNA genes, percentage of mitochondrial counts and biological sex. Linear dimensionality reduction was performed by principal component analysis (runPCA, npcs = 25). A k-nearest-neighbours graph was constructed on the basis of Euclidean distance in PCA space and refined (FindNeighbors, npcs = 25), and then the nuclei were clustered using the Louvain algorithm (FindClusters, resolution = 0.8). snRNA clusters were visualized with UMAP (runUMAP, dims = 25). To reduce the granularity of clustering, a phylogenetic tree of cluster identities was generated from a distance matrix constructed in PCA space (BuildClusterTree) and visualized as a dendrogram (PlotClusterTree). DEGs between clusters in terminal nodes of the phylogenetic tree were calculated (FindMarkers, test.use = 'wilcox', P adj < 0.05), and clusters were merged if they had fewer than 10 DEGs with the following parameters: >0.5 avg_log[FC], <10% expression in negative nuclei, and >25% expression in positive nuclei. The final de novo snRNA-seq clusters are shown in Extended Data Fig. 10c.
Inhibitory neuron clusters (Slc32a1/Gad2+) from the neonatal multiome dataset were subsequently assigned to adult BNST Esr1+ cluster labels using Seurat. Adult BNST Esr1+ clusters (as defined above) were subsetted from the adult snRNA-seq object and randomly downsampled to 5,000 nuclei. Normalization, data scaling and linear dimensionality reduction were performed with the same parameters as for neonatal and adult Esr1+ inhibitory neuron clusters. Anchor cells between adult (reference) and neonatal (query) datasets were first identified using FindTransferAnchors. Reference cluster labels, as well as the corresponding UMAP structure, were subsequently transferred to the neonatal dataset using MapQuery. Prediction scores, which measure anchor consistency across the neighbourhood structure of reference and query datasets as previously described 85 , were used to quantify the confidence of label transfer from adult to neonatal nuclei. Extended Data Fig. 10d shows the prediction scores per reference cluster and time point of nuclei mapped onto adult reference cluster labels as well as the percentage of nuclei from each de novo cluster mapped onto each adult reference cluster (prediction score > 0.5). To further validate the quality of label transfer between adult and neonatal datasets, we computed DEGs between neonatal clusters post label transfer (FindMarkers, test.use = 'wilcox', P adj < 0.05, min.diff.pct = 0.1, avg_log[FC] > 0.5) and calculated their background-subtracted, average expression (AddModuleScore) in neonatal and adult BNST Esr1+ nuclei (visualized in Extended Data Fig. 10e).
To generate pseudo-bulk, normalized ATAC bigwig tracks for each snATAC cluster, we first re-processed the cellranger ARC output BAM file for each sample using SAMtools (-q 30 -f 2) and removed duplicate reads per cell barcode using picard MarkDuplicates (BARCODE_ TAG=CB REMOVE_DUPLICATES = true). Sinto (https://timoast.github. io/sinto/) was used to split ATAC alignments for each cluster into individual BAM files using cell barcodes extracted from the Seurat object. CPM-normalized bigwig files were computed for each pseudo-bulk BAM file using DeepTools bamCoverage (--binSize 1--normalizeUsing CPM).
To analyse the neonatal multiome snATAC data, we used ArchR 88 . Separate Arrow files were created for each multiome sample, and then merged into a single ArchR project. Gene activity scores per nucleus were calculated at the time of Arrow file creation (addGeneScore-Mat = TRUE). Metadata (cluster label, sex, time and QC metrics) were transferred from the previously generated Seurat object to the ArchR project by cell barcode-matching. Dimensionality reduction was performed on the snATAC data using ArchR's iterative Latent Semantic Indexing approach (addIterativeLSI). Per-nucleus imputation weights were added using MAGIC 89 in ArchR (addImputeWeights) to denoise sparse ATAC data for UMAP visualization. Cluster-aware ATAC peak calling was performed using ArchR's iterative overlap peak merging approach (addReproduciblePeaks, groupBy = 'cluster'). Following peak calling, CISBP human motif annotations were added for each peak (addPeakAnnotation), and chromVAR deviation scores (addDe-viationsMatrix) were calculated for each motif. In addition, chromVAR was used to calculate per-nucleus deviation scores for consensus BNSTp Nfix CUT&RUN peaks. To perform neuron identity regulator analysis (Extended Data Fig. 10g), the correlation between TF RNA expression and motif deviation score was calculated for all TFs in the CISBP motif database (correlateMatrices). TFs with a correlation coefficient >0.5 and a maximum TF RNA log 2 [FC] value between each cluster in the top 50% were classified as neuron identity regulators (coloured pink in Extended Data Fig. 10g).
For visualization of gene activity and CISBP motif deviation scores ( Fig. 3c and Extended Data Fig. 10g), scores were imputed (imputeMatrix), transferred to the original Seurat object by cell barcode matching, and visualized using FeaturePlot. Signac 90 was used to generate and store peak-by-cell count matrices for each sample. snATAC markers for each cluster were calculated (FindAllMarkers, test.use = 'LR', vars. to.regress = 'nCount_ATAC', min.pct = 0.1, min.diff.pct = 0.05, logfc. threshold = 0.15). Pseudo-bulk snATAC cluster CPM was computed for each marker peak using DeepTools multiBigwigSummary and visualized with pheatmap (Extended Data Fig. 10f). Motif enrichment analysis of snATAC marker peaks for each cluster was performed using Find-Motifs. The top three enriched motifs per snATAC cluster are shown in Extended Data Fig. 10f.
To link NE-regulated loci to sex DEGs at P4 and P14 ( Fig. 3e and Extended Data Fig. 11h), we computed the Pearson correlation coefficient between sex DEG expression and NE-regulated peak accessibility for each cluster (LinkPeaks, min.distance = 2,000, distance = 1,000,000, min.cells = 2% of cluster size). Sex DEG log 2 [FC] values and NE-regulated ATAC site correlation coefficients were hierarchically clustered and visualized using ComplexHeatmap 92 .
P14 snRNA-seq data processing and analysis. Raw sequencing data were processed using the Cell Ranger pipeline (v6.0.0) with the refdata-gex-mm10-2020-A reference. Default parameters were used to align reads, count unique transcripts and filter high-quality nuclei. Individual HDF5 files for each sample were loaded into Seurat. Nuclei with lower-end RNA QC metrics (<1,000 counts) were removed. Dou-bletFinder 87 was then used to remove predicted doublets from each sample (nExp = 9% of nuclei per sample). Following doublet removal, nuclei surpassing upper-end RNA QC metrics (>20,000 counts, >6,000 genes detected) were removed. After filtering, Seurat objects were merged. Gene counts were normalized and scaled, as described for the single-nucleus multiome data processing.
The P14 snRNA-seq dataset was assigned to adult BNST inhibitory cluster labels using Seurat. Adult BNST inhibitory clusters were subsetted from the adult snRNA-seq object and randomly downsampled to 10,000 nuclei. Normalization, data scaling and linear dimensionality reduction were performed with the same parameters for P14 and adult inhibitory neuron clusters. Label transfer was then performed as described for the single-nucleus multiome data processing. Extended Data Fig. 12b shows the prediction scores of P14 nuclei mapped onto adult reference cluster labels. To validate the quality of label transfer between adult and P14 datasets, we computed DEGs between P14 clusters post label transfer, as described above, and calculated their background-subtracted, average expression (AddModuleScore) in P14 and adult BNST inhibitory clusters (shown in Extended Data Fig. 12c). Sex DEGs between control females and and control or conditional ERα KO males were calculated for each P14 cluster, as described above for the multiome analysis. Cluster abundance for each group was computed and is plotted in Extended Data Fig. 12d.
Neonatal bulk nuclear RNA-seq data processing and analysis. Reads were trimmed to remove Illumina adapters and low-quality basecalls (cutadapt -q 30), and then mapped to the mm10 reference genome using STAR. Technical duplicate reads (identical start and end positions with the same strand orientation and identical molecular identifiers) were removed using the nudup.py python package (https://github.com/ tecangenomics/nudup). The number of reads mapping to each gene (including introns) on each strand (-s 1) was calculated with feature-Counts 82 , using the mm10.refGene.gtf file. Differential gene expression analysis was performed using DESeq2 (design = ~ treatment) after prefiltering genes by expression (rowMeans ≥ 5).

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this paper.